What is this?

E11.5 (or E11?) mouse RPE

Run on 10x v3 platform

This is a trial run to determine whether Elan can process mouse tissue at this very early time point.

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.6     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.4     ✓ stringr 1.4.0
## ✓ readr   2.0.2     ✓ forcats 0.5.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(Seurat)

## Registered S3 method overwritten by 'spatstat.geom':
##   method     from
##   print.boxx cli

## Attaching SeuratObject

seurat.data  <- Read10X(data.dir = '~/data/sc_mouse_OFC_e11/aggr/outs/count/filtered_feature_bc_matrix')
seurat <- CreateSeuratObject(counts = seurat.data, project = "e11_mouse_eye_cup", min.cells = 3, min.features = 200)
seurat[["percent.mt"]] <- PercentageFeatureSet(seurat, pattern = "^mt-")

seurat@meta.data$replicate <- str_extract(colnames(seurat), '\\d')
seurat@meta.data$orig.ident <- paste0(seurat@meta.data$orig.ident, '__', seurat@meta.data$replicate)
Idents(seurat) <- seurat@meta.data$orig.ident

Basic stats per replicate

Many many more cells in rep4.

seurat@meta.data %>% group_by(orig.ident) %>% summarise(Count = n()) %>% 
  ggplot(aes(x=orig.ident,y=Count, label = Count)) + geom_bar(stat = 'identity') +
  geom_text(position =position_nudge(y=200)) +
  coord_flip() +
  cowplot::theme_cowplot()

Feature stats

But …. everything looks fairly the same between all four replicates…. so let’s just keep it all for now? There’s a larger size range in rep 4 of number of unique genes (nFeature_RNA) captured for some reason.

VlnPlot(seurat, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, group.by = 'orig.ident')

# Filtering Remove cells with percent.mt > 10%

seurat <- subset(seurat, subset = nFeature_RNA > 200  & percent.mt < 10)

seurat@meta.data %>% group_by(orig.ident) %>% summarise(Count = n()) %>% 
  ggplot(aes(x=orig.ident,y=Count, label = Count)) + geom_bar(stat = 'identity') +
  geom_text(position =position_nudge(y=200)) +
  coord_flip() +
  cowplot::theme_cowplot()

# PCA

seurat <- NormalizeData(seurat, normalization.method = "LogNormalize", scale.factor = 10000)



seurat <- FindVariableFeatures(seurat, selection.method = "vst", nfeatures = 2000)

all.genes <- rownames(seurat)
seurat <- ScaleData(seurat, features = all.genes)

## Centering and scaling data matrix

seurat <- RunPCA(seurat, features = VariableFeatures(object = seurat), npcs = 50)

## PC_ 1 
## Positive:  Dct, Gpnmb, Pmel, Tyr, Tyrp1, Bace2, Slc7a8, Krt18, Rlbp1, Trpm1 
##     Mitf, Mgll, Tspan10, Slc38a8, Oca2, Slc4a5, Gabrb3, Zdhhc2, Rab27a, S100a1 
##     Clu, Bdh2, Slc45a2, Uap1l1, Gmpr, Slc11a1, Spon1, Inmt, Ankrd44, Otx2 
## Negative:  Hmga2, Ccnd1, Pclaf, Rrm2, H19, Igfbp4, Tmsb4x, Fn1, Crabp1, Top2a 
##     Cdca8, Smc4, Vcan, Cks2, H2afx, Tubb4b, Mki67, Tubb2b, Zfhx4, Uhrf1 
##     Crym, Cenpa, Scube1, Twist1, Ccna2, Foxc1, Cenpf, Pitx2, Tpx2, Hist1h1b 
## PC_ 2 
## Positive:  Pax6, Six6, Lhx2, Sfrp2, Six3, Six3os1, Mab21l2, Ptn, Sox2, Col9a1 
##     Rax, Fign, Meis1, Syt11, Necab2, Zic1, Zic5, Car14, Hmx1, Rmst 
##     Rorb, Vsx2, Nr2f1, Dkk3, Cdon, Trpm3, Fzd5, Map7, Dapl1, Fgf15 
## Negative:  Lcp1, Fxyd5, Mef2c, Klhl6, Gmfg, Stab1, Cd34, Rasgrp3, Gngt2, Tgfb1 
##     Pecam1, Ptpn18, Vsir, Tgfbr2, Emcn, Cdh5, Cd93, Plxnd1, Hhex, Laptm5 
##     Abi3, Inpp5d, Col15a1, Smagp, Esam, Ecscr, Arhgdib, Icam1, Eng, Flt1 
## PC_ 3 
## Positive:  Ckb, Six6, Pou2f2, Sox2, Six3, Six3os1, Rax, Sfrp2, Pax6, Lhx2 
##     Zic5, Vsx2, Hmx1, Fzd5, Zic1, Emb, Meis1, Ikzf1, Rac2, Syt11 
##     Rmst, Rorb, Dkk3, C1qb, Coro1a, Dapl1, Tyrobp, Fgf15, C1qc, Cx3cr1 
## Negative:  Foxc1, Prrx1, Cald1, Clmp, Crym, Irx3, Rarb, Igf1, Col3a1, Ddr2 
##     Hmcn1, Foxc2, Parm1, Pitx2, Lum, Dkk2, Twist1, Col1a2, Lhfpl2, Rspo1 
##     Irx5, Pdgfra, Cped1, Mfap4, Col26a1, Cxcl12, Fgfr2, Foxp2, Irx1, Pdgfrb 
## PC_ 4 
## Positive:  Ikzf1, Ncf4, Blvrb, Mt1, Fth1, Acp5, Slc4a1, Gypa, Hba-a2, Alas2 
##     Rhd, Lgals1, Hemgn, Hba-a1, Ermap, Rhag, Hba-x, Hbb-y, Hbb-bt, Trim10 
##     Icam4, Hbb-bh1, Mt2, Kel, Sptb, Epb42, Prrx1, Cited4, Spta1, Car2 
## Negative:  Cdh5, Cd93, Esam, Flt1, Emcn, Tie1, Adgrf5, S100a16, Cldn5, Rasip1 
##     Col15a1, Pecam1, Sox18, Icam2, Aplnr, Ecscr, Cd34, Flt4, Vamp5, Crip2 
##     Eng, Plvap, She, Tek, Gimap6, Mmrn2, Sparcl1, Ramp2, Arhgef15, Tmem255a 
## PC_ 5 
## Positive:  Tmsb4x, Hmga2, C1qb, Coro1a, Tyrobp, Cx3cr1, C1qc, C1qa, Fcer1g, Aif1 
##     Ly86, Ctss, Adgre1, Fcgr3, Cd52, Spi1, Ptprc, Nr2f1, Trem2, Unc93b1 
##     Mpeg1, Ncf1, Fyb, Cd53, Ly6e, Ncf2, Ptpn6, Pld4, Arhgap30, Ms4a6c 
## Negative:  Slc4a1, Alas2, Gypa, Rhd, Hemgn, Rhag, Car2, Icam4, Trim10, Acp5 
##     Ermap, Hbb-bh1, Hba-a1, Hbb-bt, Epb42, Hba-a2, Hba-x, Cited4, Kel, Hbb-y 
##     Spta1, Klf1, Sptb, Pdzk1ip1, Tal1, Cldn13, Slc25a37, Hbb-bs, Slc38a5, Mgst3

print(seurat[["pca"]], dims = 1:10, nfeatures = 10)

## PC_ 1 
## Positive:  Dct, Gpnmb, Pmel, Tyr, Tyrp1, Bace2, Slc7a8, Krt18, Rlbp1, Trpm1 
## Negative:  Hmga2, Ccnd1, Pclaf, Rrm2, H19, Igfbp4, Tmsb4x, Fn1, Crabp1, Top2a 
## PC_ 2 
## Positive:  Pax6, Six6, Lhx2, Sfrp2, Six3, Six3os1, Mab21l2, Ptn, Sox2, Col9a1 
## Negative:  Lcp1, Fxyd5, Mef2c, Klhl6, Gmfg, Stab1, Cd34, Rasgrp3, Gngt2, Tgfb1 
## PC_ 3 
## Positive:  Ckb, Six6, Pou2f2, Sox2, Six3, Six3os1, Rax, Sfrp2, Pax6, Lhx2 
## Negative:  Foxc1, Prrx1, Cald1, Clmp, Crym, Irx3, Rarb, Igf1, Col3a1, Ddr2 
## PC_ 4 
## Positive:  Ikzf1, Ncf4, Blvrb, Mt1, Fth1, Acp5, Slc4a1, Gypa, Hba-a2, Alas2 
## Negative:  Cdh5, Cd93, Esam, Flt1, Emcn, Tie1, Adgrf5, S100a16, Cldn5, Rasip1 
## PC_ 5 
## Positive:  Tmsb4x, Hmga2, C1qb, Coro1a, Tyrobp, Cx3cr1, C1qc, C1qa, Fcer1g, Aif1 
## Negative:  Slc4a1, Alas2, Gypa, Rhd, Hemgn, Rhag, Car2, Icam4, Trim10, Acp5 
## PC_ 6 
## Positive:  Rbp1, Vim, Col9a1, Foxp2, Tgfb2, Ccnd1, Mfap4, Lhx2, Prrx1, Meg3 
## Negative:  Wnt4, Wnt6, Perp, Fermt1, Sfn, Esrp1, Trp63, Cldn6, Cdh1, Barx2 
## PC_ 7 
## Positive:  Ccnd2, Ung, Cp, Smoc1, Slitrk1, Ccnd1, Miat, Pax2, Cdkn1c, Nnat 
## Negative:  Cenpf, Prc1, Ube2c, Top2a, Nusap1, Cenpe, Hmmr, Cks2, Kif11, Ccnb1 
## PC_ 8 
## Positive:  Pax2, Vax1, Bmpr1b, Smoc1, Ntn1, Fzd8, Cp, Slitrk1, Adamts16, Nr2f1 
## Negative:  Rbp1, Tbx2, Dapl1, Rax, Id1, Ldha, Wfdc1, Unc5b, Jag1, Sfrp2 
## PC_ 9 
## Positive:  Cenpa, Aldh1a3, Ccnb2, Cdca3, Tubb4b, Trpm1, H2afx, Ube2c, Cks2, AY036118 
## Negative:  Kcnq1ot1, Map1b, Epha5, Vcan, Sptbn1, Ncam1, Ank2, Itm2a, Abca5, Prss23 
## PC_ 10 
## Positive:  Filip1, Adamts9, Hcn1, Akap13, Tbx20, Fam171b, Gja1, Dazl, Mab21l2, H19 
## Negative:  Col4a4, Crabp1, Gm534, Gm33489, Insm1, Ttr, Espn, Vit, Olfm1, Col8a1

DimPlot(seurat, reduction = "pca")

# Heatmaps by PC

DimHeatmap(seurat, dims = 1:30, cells = 500, balanced = TRUE)

# Jackstraw plot

Goal is to roughly figure out how many Principal Components (from PCA) to use downstream. I picked 20.

seurat <- JackStraw(seurat, num.replicate = 100, dims = 50)
seurat <- ScoreJackStraw(seurat, dims = 1:50)


JackStrawPlot(seurat, dims = 1:50)

## Warning: Removed 74068 rows containing missing values (geom_point).

ElbowPlot(seurat, ndims = 50)

# Clustering and get markers for each cluster

seurat <- FindNeighbors(seurat, dims = 1:20)

## Computing nearest neighbor graph

## Computing SNN

seurat <- FindClusters(seurat, resolution = 0.5)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 4034
## Number of edges: 119026
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8926
## Number of communities: 12
## Elapsed time: 0 seconds

seurat <- RunUMAP(seurat, dims = 1:20)

## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session

## 13:47:23 UMAP embedding parameters a = 0.9922 b = 1.112

## 13:47:23 Read 4034 rows and found 20 numeric columns

## 13:47:23 Using Annoy for neighbor search, n_neighbors = 30

## 13:47:23 Building Annoy index with metric = cosine, n_trees = 50

## 0%   10   20   30   40   50   60   70   80   90   100%

## [----|----|----|----|----|----|----|----|----|----|

## **************************************************|
## 13:47:24 Writing NN index file to temp file /var/folders/s4/y5f1tt296dj8088gvczcx11d4lrnr7/T//RtmpXS9RWv/file15416691101d1
## 13:47:24 Searching Annoy index using 1 thread, search_k = 3000
## 13:47:25 Annoy recall = 100%
## 13:47:25 Commencing smooth kNN distance calibration using 1 thread
## 13:47:26 Initializing from normalized Laplacian + noise
## 13:47:26 Commencing optimization for 500 epochs, with 154014 positive edges
## 13:47:30 Optimization finished

DimPlot(seurat, reduction = "umap")

seurat.markers <- list()
for (i in seurat@meta.data$seurat_clusters %>% unique()){
  seurat.markers[[i]] <- FindMarkers(seurat, ident.1 = i, only.pos = TRUE, logfc.threshold = 0.5)
}

Regress out G2M <-> S phase diff

As the previous analysis essentially had clusters by cell phase, which is of less interest

s.genes <- cc.genes$s.genes
g2m.genes <- cc.genes$g2m.genes


m.s.genes <-s.genes %>% str_to_title()
m.g2m.genes <- g2m.genes %>% str_to_title()
seurat <- CellCycleScoring(seurat, s.features = m.s.genes, g2m.features = m.g2m.genes)

## Warning: The following features are not present in the object: Mlf1ip, not
## searching for symbol synonyms

## Warning: The following features are not present in the object: Fam64a, Hn1, not
## searching for symbol synonyms

DimPlot(seurat, group.by = 'Phase')

seurat$CC.Difference <- seurat$S.Score - seurat$G2M.Score
options(future.globals.maxSize = 8000 * 1024^2)
seuratCC <- ScaleData(seurat, vars.to.regress = "CC.Difference", features = rownames(seurat))

## Regressing out CC.Difference

## Centering and scaling data matrix

seuratCC <- RunPCA(seuratCC, features = VariableFeatures(object = seuratCC), npcs = 50)

## PC_ 1 
## Positive:  Dct, Pmel, Gpnmb, Tyr, Tyrp1, Bace2, Slc7a8, Krt18, Rlbp1, Trpm1 
##     Mgll, Mitf, Tspan10, Slc38a8, Oca2, Slc4a5, Gabrb3, Zdhhc2, Rab27a, S100a1 
##     Clu, Bdh2, Uap1l1, Slc45a2, Gmpr, Slc11a1, Spon1, Inmt, Ankrd44, Ugcg 
## Negative:  Hmga2, Ccnd1, Rrm2, Pclaf, Cks2, H19, Cdca8, Igfbp4, Tmsb4x, Uhrf1 
##     Crabp1, Top2a, Fn1, H2afx, Cenpf, Smc4, Vcan, Ube2c, Tubb4b, Mki67 
##     Cenpa, Tpx2, Pcna, Tubb2b, Zfhx4, Dtl, Ccna2, Crym, Cdca3, Scube1 
## PC_ 2 
## Positive:  Pax6, Lhx2, Six6, Sfrp2, Six3, Mab21l2, Six3os1, Ptn, Sox2, Col9a1 
##     Rax, Fign, Meis1, Necab2, Syt11, Car14, Zic1, Zic5, Hmx1, Rmst 
##     Rorb, Nr2f1, Vsx2, Dkk3, Trpm3, Cdon, Map7, Fzd5, Dapl1, Fgf15 
## Negative:  Lcp1, Fxyd5, Mef2c, Klhl6, Gmfg, Stab1, Cd34, Rasgrp3, Gngt2, Tgfb1 
##     Pecam1, Ptpn18, Vsir, Tgfbr2, Emcn, Cdh5, Cd93, Plxnd1, Hhex, Laptm5 
##     Abi3, Inpp5d, Col15a1, Esam, Smagp, Ecscr, Arhgdib, Icam1, Eng, Flt1 
## PC_ 3 
## Positive:  Ckb, Six6, Pou2f2, Six3, Sox2, Six3os1, Rax, Sfrp2, Pax6, Lhx2 
##     Zic5, Ikzf1, Vsx2, Hmx1, Fzd5, Zic1, Emb, Meis1, Rac2, Rmst 
##     Syt11, C1qb, Rorb, Coro1a, Dkk3, Tyrobp, C1qc, Fcer1g, Cx3cr1, Dapl1 
## Negative:  Foxc1, Prrx1, Cald1, Clmp, Crym, Irx3, Rarb, Igf1, Col3a1, Hmcn1 
##     Ddr2, Foxc2, Parm1, Lum, Pitx2, Dkk2, Col1a2, Twist1, Lhfpl2, Rspo1 
##     Irx5, Cped1, Mfap4, Pdgfra, Col26a1, Cxcl12, Foxp2, Fgfr2, Irx1, Pdgfrb 
## PC_ 4 
## Positive:  Ncf4, Blvrb, Ikzf1, Mt1, Fth1, Acp5, Slc4a1, Gypa, Hba-a2, Alas2 
##     Rhd, Hemgn, Hba-a1, Hba-x, Hbb-y, Rhag, Ermap, Hbb-bt, Trim10, Icam4 
##     Hbb-bh1, Lgals1, Kel, Sptb, Mt2, Epb42, Spta1, Cited4, Car2, Klf1 
## Negative:  Cdh5, Cd93, Esam, Flt1, Emcn, Tie1, Adgrf5, S100a16, Cldn5, Rasip1 
##     Col15a1, Pecam1, Icam2, Sox18, Aplnr, Cd34, Ecscr, Flt4, Crip2, Eng 
##     Vamp5, She, Plvap, Tek, Mest, Gimap6, Mmrn2, Sparcl1, Arhgef15, Tmem255a 
## PC_ 5 
## Positive:  Tmsb4x, C1qb, Tyrobp, Coro1a, Cx3cr1, C1qc, C1qa, Hmga2, Fcer1g, Aif1 
##     Ly86, Ctss, Adgre1, Fcgr3, Cd52, Spi1, Ptprc, Trem2, Unc93b1, Mpeg1 
##     Fyb, Ncf1, Cd53, Ncf2, Pld4, Ptpn6, Arhgap30, Csf1r, Ms4a6c, C3ar1 
## Negative:  Slc4a1, Alas2, Gypa, Rhd, Hemgn, Rhag, Car2, Icam4, Trim10, Hbb-bh1 
##     Acp5, Hba-a1, Ermap, Hbb-bt, Hba-x, Hba-a2, Epb42, Cited4, Kel, Hbb-y 
##     Spta1, Klf1, Sptb, Pdzk1ip1, Tal1, Cldn13, Slc25a37, Hbb-bs, Slc38a5, Mgst3

DimPlot(seuratCC, reduction = "pca")

seuratCC <- FindNeighbors(seuratCC, dims = 1:20)

## Computing nearest neighbor graph

## Computing SNN

seuratCC <- FindClusters(seuratCC, resolution = 0.5)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 4034
## Number of edges: 122396
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8868
## Number of communities: 14
## Elapsed time: 0 seconds

seuratCC <- RunUMAP(seuratCC, dims = 1:20)

## 14:45:10 UMAP embedding parameters a = 0.9922 b = 1.112

## 14:45:10 Read 4034 rows and found 20 numeric columns

## 14:45:10 Using Annoy for neighbor search, n_neighbors = 30

## 14:45:10 Building Annoy index with metric = cosine, n_trees = 50

## 0%   10   20   30   40   50   60   70   80   90   100%

## [----|----|----|----|----|----|----|----|----|----|

## **************************************************|
## 14:45:11 Writing NN index file to temp file /var/folders/s4/y5f1tt296dj8088gvczcx11d4lrnr7/T//RtmpXS9RWv/file1541665f408cb
## 14:45:11 Searching Annoy index using 1 thread, search_k = 3000
## 14:45:11 Annoy recall = 100%
## 14:45:12 Commencing smooth kNN distance calibration using 1 thread
## 14:45:12 Initializing from normalized Laplacian + noise
## 14:45:12 Commencing optimization for 500 epochs, with 155152 positive edges
## 14:45:17 Optimization finished

DimPlot(seuratCC, reduction = "umap")

DimPlot(seuratCC, reduction = "umap", group.by = 'Phase')

seurat.markersCC <- list()
for (i in seuratCC@meta.data$seurat_clusters %>% unique()){
  seurat.markersCC[[i]] <- FindMarkers(seuratCC, ident.1 = i, only.pos = TRUE, logfc.threshold = 0.5)
}

Output files

save(seurat.markers, seurat, seurat.markersCC, seuratCC, file = '../data/00_process_output.Rdata')

Session

sessionInfo()